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In light filamentation induced by axicon-generated, powerful Bessel beams, the spatial propagation 
dynamics in the nonlinear medium determines the geometry of the filament channel and hence its 
potential applications. We show that the observed steady and unsteady Bessel beam propagation 
regimes can be understood in a unihed way from the existence of an attractor and its stability 
properties. The attractor is identified as the nonlinear unbalanced Bessel beam (NL-UBB) whose 
inward Hankel beam amplitude equals the amplitude of the linear Bessel beam that the axicon would 
generate in linear propagation. A simple analytical formula that determines de NL-UBB attractor 
is given. Steady or unsteady propagation depends on whether the attracting NL-UBB has a small, 
exponentially growing, unstable mode. In case of unsteady propagation, periodic, quasi-periodic or 
chaotic dynamics after the axicon reproduces similar dynamics after the development of the small 
unstable mode into the large perturbation regime. 


I. INTRODUCTION 

The nonlinear propagation of very intense pulsed 
Bessel beams (BBs) has attracted a lot of attention in 
recent years, specially because of the ability of BBs of 
creating filamentary ionized channels that may be longer 
and more spatially controllable n that the filaments 
created by focusing standard, Gaussian-like light pulses 
Si- The versatility of Bessel beams for filamentation 
has been dramatically demonstrated very recently with 
the generation of tubular plasma channels when the seed¬ 
ing Bessel beam carries an optical vortex ii- These 
achievements have opened new perspectives in ultrafast 
laser material processing in transparent dielectrics, such 
as waveguide writing and micro- or nanomachining [^, Q , 
or in long-range filamentation in gases with application, 
for instance, in microwave guiding by filaments in the 
atmosphere Q. 

The first studies on BB propagation in nonlinear media 
date from the beginning of the past decade [1, [l^ ■ Non¬ 
linear Bessel beams as stationary (non-diffracting) solu¬ 
tions to the nonlinear Schrodinger equation (NLSE) with 
Kerr nonlinearity were first introduced in 11 ill . Nonlin¬ 
ear unbalanced Bessel beams (NL-UBBs) [l^ were later 
found as stationary solutions of the NLSE in media with 
Kerr nonlinearity and nonlinear losses (NLLs), just the 
two key nonlinearities determining the spatial dynamics 
in BB filamentation. These NL-UBBs have indeed been 
proven to play a prominent role in the filamentation with 
axicon focused BBs [I|, [H, [l3| , acting as attractors of the 
dynamics. Matter waves of this kind can also exist in 
Bose-Einstein condensates [l^ . More recently NL-UBBs 
carrying vortices have been also described [l6| , and have 
similarly found to act as attractors in the filamentation 
seeded by axicon focused vortex BBs 


In the experimental and numerical studies on nonlin¬ 
ear BB propagation, two different initial conditions for 
the light entering the medium are usually considered. In 
a first arrangement, BBs are launched into the medium 
when they are already formed e. g., 

an ideal BB at any transversal plane, or an apodized 
BB at the focus of an axicon. Except if NLLs domi¬ 
nate initially the dynamics [l6l - [T8l| . Kerr nonlinearity in¬ 
duces in this case large temporal and spatial instabilities 
l- d 1911 ■ In most of filamentation experiments with BBs 
Ini ■ Il3l [T3| and related numerical studies [2^ , the radi¬ 
ation exiting from the BB generator enters the nonlinear 
medium prior to the formation of the BB, in a state of 
widespread energy at low intensity levels, so that the lin¬ 
ear BB is never formed. With an axicon, for example, the 
medium is placed in contact with it, or simply fills the 
space surrounding it, as in filamentation in gases. This 
“soft” input condition has been proven useful to prevent 
the onset of large temporal instabilities in the nonlinear 
medium 0. With this arrangement, two different Bessel 
beam propagation regimes have been observed [10. 
In a steady Bessel propagation regime, the input radia¬ 
tion undergoes a transformation into a quasi-stationary 
state within the Bessel zone that has been identified as 
a NL-UBB. In a unsteady regime, the light intensity and 
fluence feature periodic, quasi-periodic, even disordered 
spikes in the Bessel zone (and azimuthal breaking in the 
case of vortex BBs i). This regime has been associated 
with sufficiently small cone angles and relatively low in¬ 
put powers so that self-focusing is the dominant nonlin¬ 
earity. 

In this paper we aim at providing a unified understand¬ 
ing of these two regimes of BB propagation under soft 
input conditions. We show that these two regimes are 
different manifestations of the same underlying dynam¬ 
ics. Either steady or unsteady, the spatial dynamics is 
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dominated by the existence of an attractor in the form of 
a specific NL-UBB. We identify the attracting NL-UBB 
and derive an approximate analytical expression specify¬ 
ing it in terms of the properties of the nonlinear medium 
and the light beam illuminating the axicon. However, an 
attractor is not necessarily a stable attractor; its insta¬ 
bility may lead to a richer dynamics around it, includ¬ 
ing periodic, quasi-periodic and chaotic behavior. The 
unsteady or steady regimes are seen to be determined 
by the existence or not of a small, exponentially grow¬ 
ing unstable mode of the attracting NL-UBB. In the un¬ 
steady regime, the unstable dynamics in the Bessel zone 
of the axicon is triggered by the unstable mode and is 
seen to reproduce its characteristic oscillation frequency, 
its development into large periodic or quasi-periodic an- 
harmonic oscillations, or its development into chaotic 
oscillations, depending on the gain of the small unsta¬ 
ble mode. Although the unsteady Bessel filamentation 
regime has been previously s ugg ested to be associated 
with NL-UBB instability [l^ Il4l | . it is only the identifi¬ 
cation of the attracting NL-UBB that have allowed us to 
analyze its stability properties, and hence to verify that 
hypothesis, putting it in quantitative terms. 

For simplicity we focus on BBs generated by axicons in 
most of the numerical simulations, but the same results 
are seen to hold for other soft input conditions that would 
generated BBs in linear propagation. We illustrate the 
results in air at 800 nm, in which case the characteristic 
angles separating the different regimes are quite small, 
but we have verified that the same results hold at larger 
angles (but still paraxial) in condensed media. 


II. NONLINEAR UNBALANCED BESSEL 
BEAMS 

We consider diffraction, Kerr nonlinearity and NLLs 
as the key effects determining the propagation of the 
light beam coming from the BB generator. In the parax¬ 
ial approximation, the envelope A of the light beam 
E = Aexp(—iujt + ikz) of frequency w and propagation 
constant k, is then suitably described by the NLSE 

= + ( 1 ) 

where n, n 2 and are, respectively, the linear and 

nonlinear refractive indexes and the M-photon absorp¬ 
tion coefficient. For the initial conditions of interest, and 
according to [3 j temporal effects are assumed to play a 
secondary role. 


In order to properly understand the propagation, it is im¬ 
portant to review the properties of NL-UBBs, stressing 
their asymptotic properties. NL-UBBs were introduced 
in [T^ as non-diffracting and non-attenuating solutions 
of (HI) of the form A = a(r) exp[*(/)(r)] exp(i(52:), where 



FIG. 1. (a) In air at 800 nm (n ~ 1, n 2 = 3.2 x 10”^^ 

cm^/W, M = 8, = 1-8 X 10“®'^ cm^®/W^), region in pa¬ 

rameters space {Io,0) of existence of NL-UBBs. The dashed 
black curve is the approximate analytical curve (lo.max,^) 
given in the text. To the right of the gray dotted curve, 
NL-UBB are said to be NLL-dominated in the sense that 
Lmpa < Ldif and Lmpa < I/Kerr, where the characteristic 
lengths are Lmpa = 2 /for multiphoton absortion, 
Ldif = 1/|<5| (Rayleigh range of central lobe of Bessel function) 
for diffraction, and I/Kerr = n/kn 2 lo for Kerr nonlinearity. 
Insets 1 and 2: Amplitude radial profiles of NL-UBBs with 
6 = 0.15 deg, and lo = 12 and 28 TW/cm^ (solid curves), 
and of linear BBs of the same intensity (dashed curves). The 
squares 1 and 2 indicate the location of those NL-UBBs in pa¬ 
rameters space, (b) Values of |6in P and |6out|^ for NL-UBBs 
with 9 = 0.15 and 1 degrees, extracted from the numerically 
evaluated radial profiles (gray dotted curves), and approxi¬ 
mately evaluated from Eq. ([§]) (solid curves). 


S = —k6^j2 is a negative axial wave vector shift corre¬ 
sponding (in the paraxial approximation) to a cone angle 
9. The real amplitude a{r) > 0 and phase (j){r) satisfy 

a"+ ^ + = 0, ( 2 ) 

r n 

—Fr =- 2TTr(j)'a^ = f3^^'^2TT [ drra^^ = Nr , (3) 

k Jo 

(prime signs stand for d/dr) with boundary conditions 
a(0) = ^/Io, a'(0) = 0, </>'(0) = 0, where /q is the 
peak intensity of the NL-UBB. Equation ([3]) is the re¬ 
filling condition for stationarity with nonlinear absorp¬ 
tion, stating that the nonlinear power losses W in each 
circle of radius r are compensated by an inward radial 
flux —Fr through its circumference. For each cone an¬ 
gle 0, NL-UBB exist up to a maximum value /o,max 
of the peak intensity whose value depends on the op¬ 
tical properties of the medium at w. Figure [Ija) shows 
the region of existence in the parameters space {Io,9) 
of NL-UBBs in air at 800 nm, and the insets two typ¬ 
ical radial amplitude profiles. An approximate for¬ 
mula relating Jo,max to 9 is 9^ = I^^^Jko - 

2n2Jo,max/no [dashed curve in Fig. (Ha)], where 
ctm = 0.542,0.381,0.313,0.295,0.244,0.223... for M = 
3,4,... [l^. At large radius r, NL-UBBs behave asymp¬ 
totically as BBs, but with unbalanced amplitudes of its 
outward and inward Hankel components: 


~ 1 [5outiJ,$'^(fc0r) + 6miJ^")(fc0r) 


(4) 
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with |6in| > |6out|, and only bi^ = bout = \/Tb for a lin¬ 
ear BB of intensity Ib- The amplitudes |&in| and |&out| 
can be easily extracted from the radial intensity profile, 
known numerically or experimentally: Using the asymp¬ 
totic forms of Hankel functions for large argument one 
gets, from Eq. dU, ~ (a^/r)[l -I- Ccos{2k9r -\- tp)], 
with C = 2|6out||^in|/(|^outP + l^inp) and = (|6outP + 
|^inP)/(27r), meaning that the asymptotic radial intensity 
proHle consists of oscillations of contrast C < 1 about an 
average value a^lr. From these features, the Hankel am¬ 
plitudes are obtained to be |5in_outP = 7ra^(l±Vl — U^)- 
Examples of their values extracted from the numerical in¬ 
tensity profiles are depicted in Fig. [TJb) for increasing 
NL-UBB intensities Iq at two cone angles. 


III. ON THE DYNAMICS OF REAL BESSEL 
BEAMS IN NONLINEAR MEDIA 

In Ref. [T^, it was demonstrated experimentally that 
a BB ^(r, 0) = ^/IsJoikOr) launched in a nonlinear 
medium in a regime where NLLs are significant trans¬ 
forms spontaneously into a NL-UBB that preserves the 
cone angle. In the ideal case of BBs carrying infinite 
power (and for the broader class of vortex NL-UBBs) the 
specific attracting NL-UBB has been identified as that 
whose inward Hankel amplitude equal the amplitude of 
the launched BB, that is, |6in \ = VT^M- Since for the 
input BB bin = bout = \/Tb, the amplitude of the inward 
Hankel component can be said to be a preserved quantity 
in the nonlinear dynamics. 

In actual settings the input power is finite, and the 
medium is placed close to or in contact with an axicon 
(or other BB generators), or simply fills the space sur¬ 
rounding the axicon, so that the BB is not formed when 
the radiation enters the medium and propagates linearly 
initially. With an axicon, for instance, the field entering 
the medium at z = 0 is usually modelled by 

A(r, 0) = exp(—exp(—(5) 

where w and Iq are the width and the intensity of the 
Gaussian beam illuminating the axicon. In linear prop¬ 
agation, this would produce an apodized BB of inten¬ 
sity Ib = TTkw6Ic/y/e at a distance zb = wj20, which 
is one-half the length w/9 of the so-called Bessel zone. 
Under these soft input conditions, the unsteady Bessel 
filamentation regime, associated with small cone angles 
and relatively low intensities, results in periodic or quasi- 
periodic field oscillations [l|, [l^ , though it can also result 
in chaos (see below). The steady filamentation regime, 
associated with large cone angles or higher intensities, 
has been explained in terms of the formation of a NL- 
UBB [llildl. 

As pointed out in the introduction, these regimes are 
shown here to be different manifestations of the same un¬ 
derlying dynamics. Either steady or unsteady the spa¬ 
tial dynamics in BB filamentation is dominated by the 


existence of an attractor in the form of a specific NL- 
UBB, unsteady or steady regimes being determined by 
the existence or not of a small, exponentially growing 
unstable mode of the attracting NL-UBB. In Section ITVl 
we identify the attractor and obtain approximate analyt¬ 
ical formulas determining it. In Section]^ we perform a 
linearized stability analysis of ideal NL-UBB and note a 
biunivocal relation between NLUBB stability/insability 
under small perturbations and steady/unsteady propa¬ 
gation after the real BB generator. The unsteady Bessel 
regime appears then to be triggered by the existence of 
a small unstable mode in the attracting NL-UBB. De¬ 
pending on the gain, signatures of this mode, or of its de¬ 
velopment into large periodic, quasi-periodic or chaotic 
perturbation regimes are indeed observed in the Bessel 
zone of the axicon. 


IV. THE ATTRACTING NONLINEAR 
UNBALANCED BESSEL BEAM 

We identify the attracting NL-UBB as that whose in¬ 
ward Hankel amplitude coincides with the amplitude of 
the BB that the BB generator would create in linear prop¬ 
agation, i. e., the NL-UBB with |&in| = ^/TB■ Since 
bin = bout = for BBs, the amplitude of the inward 
Hankel component is not affected by nonlinearities, and 
in this sense can be said to be conserved. This conclu¬ 
sion is extracted from extensive numerical simulations, 
of which only a few examples are shown. Conceptually, 
it is not difficult to understand that the inward Hankel 
component created by the axicon, even if of finite power, 
and supplying power conically inwards is not affected by 
nonlinear absorption at the beam center in the Bessel 
zone. 


Figure [2] illustrates this law for an axicon imprint¬ 
ing a cone angle 9 = 0.15 deg illuminated by Gaussian 
beams of width w = 1.5 cm and intensities Iq = 0.0666 
and Iq = 0.0174 TW/cm^ and propagating in air at 
A = 800 nm. In linear propagation these Gaussian 
beams would create linear BBs [dashed curves in Figs. 
UKa) and (c)] of intensities Ib = 39.17 TW/cm^ and 
Ib = 10.22 TW/cm^, respectively [dashed horizontal 
lines Ib in Figs. [2Ka) and (c)]. In the steady regime of 
Fig. dKa), the on-axis intensity [solid curve in Fig. dKa)] 
in the Bessel zone stabilizes in the numerically evaluated 
intensity Iq = 28 TW/cm^ (solid horizontal line Iq) of 
the NL-UBB having |6in| = VTb = 39.17 [see Fig. HKb)]. 
The whole beam transforms in fact into the attracting 
NL-UBB, as seen in Fig. [2][b) showing intensity profiles 
at increasing distances up to the center of the Bessel zone 
Zb = 286 cm. In the unsteady regime of Fig. [2Kc), the on- 
axis intensity [solid curve in Fig. dKc)] also approaches. 
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FIG. 2. (a) For a Gaussian beam of with w = 1.5 cm and 
peak intensity Iq = 0.0666 TW/cm^ illuminating an axicon 
that forms a BB of cone angle 0 = 0.150 deg and peak inten¬ 
sity Ib = 39.17 TW/cm^ at the center zb — 286 cm of the 
Bessel zone in linear propagation (dashed curve and horizon¬ 
tal dashed line), on-axis intensity in nonlinear propagation 
in air at 800 nm (solid curve), and numerically evaluated in¬ 
tensity Jo = 28 TW/cm2 of the NL-UBB with |bin|^ = Ib 
(horizontal solid line), (b) Radial intensity profile at increas¬ 
ing propagation distances up to the center of the Bessel zone 
(solid curves), and radial intensity profile of the attracting 
NL-UBB with |6in|^ = Ib (gray dashed curve), (c) and (d) 
The same as in (a) and (b) but with la = 0.0174 TW/cm^ for 
the input Gaussian beam, Ib = 10.22 TW/cm^ for the linear 
BB, and h = 12 TW/cm^ for the NL-UBB with \binf = Ib- 


but now oscillates about the numerically evaluated in¬ 
tensity Iq = 12 TW/cm^ (solid horizontal line Iq) of the 
NL-UBB with |6in| = = 10.22 TW/cm^ [see Fig. 

HTb)], and the same happens to the whole radial profile 
at increasing distances. Oscillations may be much more 
pronounced and disordered, but clear signatures of the 
attracting NL-UBB, its dominant small unstable mode, 
and its development into a large perturbation regime, are 
always observable, as shown below (see Fig. El). 

The law |tiin| = ^/Tb holds for other finite-power ver¬ 
sions of BBs, as the Bessel-Gauss beam 

A{r,z) = 

xexp(^--^j, (6) 
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FIG. 3. The same as in Fig. [5] (a) and (c) but the beam 
entering into the medium is the Bessel-Gauss beam of Eq. 
® with V = 2 cm at 2 = —929 cm from the waist at 2 ; = 0. 
In (a), Jfl = 39.17 TW/cm^ leading to Iq = 28 TW/cm^, and 
in (b), Ib = 10.22 TW/cm^ leading to 7o = 12 TW/cm^ 


producing at z = 0 the Gaussian-apodized BB 

\/7BJo{k0r)e~^ . For a soft input into the medium, 

the entrance plane is at z = Zin 0 such that the in¬ 
tensity is low enough for nonlinear effects to be initially 
negligible. For v = 2 cm and Zin = —929 cm, propagation 
in air at the same wave length and the same linear BB 
peak intensities Ib = 39.17 and Ib = 10.22 TW/cm^ as 
in Fig. |2l the attracting NL-UBBs are seen in Fig. |3]to 
have the same intensities Iq = 28 and Iq = 12 TW/cm^ 
as with the axicon. 

Thus, given the intensity Ib that a BB generator would 
create, it is possible to foresee the attracting NL-UBB. In 
practice, this requires to extract the values of |6inP from 
the numerical radial profiles of NL-UBBs of different in¬ 
tensities Iq with the given cone angle in the particular 
medium, as explained above [dotted curve |6inP in Fig. 
EKb)] and in Ref. [l^, and to pick up the particular NL- 
UBB with |6inP = Ib- This long numerical procedure 
would be greatly simplified if we had analytical expres¬ 
sions for |6inP as functions of the NL-UBB parameters 
{9, Iq) and the optical properties of the medium. 

An approximate expression can be obtained as fol¬ 
lows. We first note that Eg . O implies —Foe = 
(l^inP — |^outP)/fc = Noo [U, meaning that the un¬ 
balance of the Hankel amplitudes sets a net constant 
inward radial power flux coming from a reservoir at 
large radial distances to refill the total NLL Noo dur¬ 
ing the propagation. At the same time, most of NLL 
take place in the beam center, where the NL-UBB can 
be approached, if the NL-UBBs is not well within the 
NLL-dominated region (see caption of Fig. [T]), by A ~ 
y/jQ Jni^/k'^O'^ + 2fcA:NT. r) exp(f Jz), with /cnl = kn 2 lQ/n 
jl^ . Evaluation of Nao with this profile yields the fol¬ 
lowing approximate expression for the NLL of a NL- 
UBB, and hence an approximate relation between |6in| 
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and |6out|: 


kNoo = \hnf - \boutf ^ 


p{M)iM 


fc6l2(l + 2n2/o/n6l2) 


v(^) 


(7) 


where 7 ^^^) = 27r Jq^{ x)xdx is a number. 

Second, numerical evaluation of |6in,out| reveals that, 
except for NLL-dominated NL-UUBs (see caption of Fig. 
[T|), their average value (|5out| + |^in|)/2 can be approxi¬ 
mated by the value |6in| = |&out| = I^Kerrl of the asymp¬ 
totic form of the solutions of Eqs. , © and m in the 
absorption-less case (i. e., with = 0). Without ab¬ 
sorption, (^(r) = 0 and the amplitude of nonlinear BBs 
behaves as a{r) ~ ^[bKeriHo^\k0r) + at 

large r. The scaling p = kOr and d = ajy/lQ in Eqs. @, 
leads to the one-parameter problem d" -I- a! jp + pa? = 
0 , with initial conditions d(0) = 1, d'(0) = 0, and 
where p = 2n2lQln6'^. Erom the numerical solution 
of this problem with different values of p, we find the 
value of I^Kerrl of the scaled asymptotic form d{p) ~ 
i[6Kerr77o^^(p) -f (p)] a function of p, which 

is found to fit accurately to the function |6Kerr| = f{p) = 
(1 -I- cp)/{l + dp) with c = 0.63 and d = 0.76. Com¬ 
ing back to real variables, I^Kerrl = + cp)/{l -I- dp) 

provides semi-analytical solution to the asymptotic form 
of nonlinear BBs in transparent media. Finally, since 
(|^out| + |^in|)/2 — I^Kerrl in the nonlinearly lossy medium, 
we obtain 


|dout| + I din I — 2 


1 + c{2n2h / nO'^) 
1 -I- d{2n2lo/n9‘^) 


'/h ■ 


( 8 ) 


The two relations © and m lead to the approximate 
formulas 


Idin.outI - f{p)\/%±P^^'^ 




Ak0^1+p)fip)^' 


(9) 


for the amplitudes of the inward and outward Hankel 
components of NL-UBBs as functions of their cone angle 
0 and peak intensity Iq and the medium properties [see 
Fig. mb), solid curves]. 

If we now set jdinj = V^Tb in Eq. ®, we obtain the 
approximate equation 


- /(p)\A) + 7 ^^^ 




4fc6»2(l +p)f{p)^ 


( 10 ) 


{p = 277 , 2/0 /nd"^) relating the intensity Iq of the attracting 
NL-UBB to the cone angle 0 and the intensity Ib of the 
linear BB that the Bessel-beam generator would create. 
As an example, the values of Iq provided by Eq. m 
for NL-UBBs in air at 800 nm, two cone angles 0 and 
increasing Ib are plotted in Eig. [d] and are seen to match 
quite accurately the numerically obtained values for NL- 
UBBs of intensities Iq below the NLL-dominated case 
(horizontal dotted lines). Equation (fTUl) is seen to give a 
reasonably good estimate of Iq even at huge intensities of 
Ib = 100 TW/cm^ well within the NLL-dominated case. 



FIG. 4. In air at 800 nm and for the indicated cone angles, 
intensity Iq of the attracting NL-UBB as a function of the in¬ 
tensity Ib of the linear BB that the BB generator would cre¬ 
ate in linear propagation, numerically evaluated (dotted gray 
curves), and obtained from Eq. (1101) (solid curves). Above 
the horizontal dotted lines NL-UBBs are NLL-dominated. 





FIG. 5. (a) Gain —Iiuk and (b) oscillation frequency Rere of 
the most unstable mode mode of NL-UBBs in air at 800 nm 
as functions of their cone angle for a few values of their peak 
intensity Iq. ( c ) and (d) The same as in (a) and (b) but as 
functions of the peak intensity Iq for a few values of the cone 
angle 6. 


V. STEADY AND UNSTEADY REGIMES 
VERSUS NL-UBB STABILITY 


Once the attractor is specified, our numerical simu¬ 
lations indicate that there is a bi-univocal relation be¬ 
tween steady/unsteady propagation regime and stabil¬ 
ity/instability of the attracting NL-UBB against small 
radial perturbations. 
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FIG. 6. In air at 800 nm, on-axis intensity I versus propagation distance z of perturbed, ideal NL-UBBs of peak intensities 
intensities lo = 14 TW/cm^ and cone angles (a) 6 = 0.075 deg and (b) 6 = 0.05 deg. (c) and (d) Corresponding phase spaces 
I-dl/dz in [400,3000] cm (beginning of the large perturbation regime), (e) and (f) Corresponding phase spaces for the whole 
propagation range [0,12000] cm. The dashed lines locate the attractor. 


For fundamental (vortex-less) NL-UBBs, radial insta¬ 
bility appears to be the dominant instability, since no az¬ 
imuthal breakiim has been observed in experiments and 
simulations [U, flM IT^ , particularly under soft in¬ 

put conditions (l4| . Linearized stability analysis of NL- 
UBB against radial perturbations has been performed 
numerically for typical values of NL-UBB parameters in 
Ref. [l2| and [ij, where all details of the procedure are 
explained. In short, supposing a solution to the NLSE 
(P) of the form 

A = a(r)e*'^(’')e*^^ -f e[u(r)e“^ -b , (11) 


that is, a NL-UBB plus a small (e —?> 0) mode (u, v) that 
grow exponentially in case that ImK < 0 while (possibly) 
oscillating harmonically with frequency Rck, a differen¬ 
tial eigenvalue problem is obtained for the eigenvalues 
K and eigenmodes modes {u,v), which has to be solved 
numerically. As pointed out in [l^ . the difficulty with 
this analysis for NL-UBBs compared to that for stan¬ 
dard solitons lies in the weak localization of NL-UBBs. 
Truncation of the NL-UBB in any finite radial box im¬ 
posed by the numerical procedure sets a lower bound to 
the reliable values of jlmAtj [l^. Figure [S] shows exam- 
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FIG. 7. (a) On axis intensity / (solid curve) in air at 800 nm after and axicon imprinting a cone angle 9 = 0.075 deg illuminated 
by a Gaussian beam of width w = 2.5 cm and intensity la = 0.02146 TW/cm^. In linear propagation (dashed curve) the 
intensity of the BB would be Ib ~ 10.508 TW/cm^ (horizontal dashed line) at zb = 955 cm, so that the attracting NL-UBB is 
defined by Iq = 14. TW/cm^ (horizontal solid line) and 6 = 0.075 deg. (b) The same except that 9 = 0.05 deg and Ig = 0.03113 
TW/cm^. In linear propagation the BB of intensity Ib = 10.165 TW/cm^ would be formed at zb = 1432 cm, so that the 
attracting NL-UBB is defined by /q = 14 TW/cm^ and 9 = 0.05 deg. (c) Phase space I - dl/dz about the focal region of the 
axicon in case (a), (d) Phase space in the whole propagation in case (b). The dashed lines in (c) and (d) locate the attractor. 


pies of the exponential gain —ImK and the oscillation 
frequency Rew of the most unstable mode of NL-UBBs 
in air at 800 nm as functions of the cone angle and fixed 
values of the peak intensity [Figs. jSKa) and (b)], and as 
functions of intensity and fixed values of the cone angle 
[Figs. [SKc) and (d)]. As noted in [Ij, NL-UBBs tend to 
stabilize as the cone angle increases. According to our 
analysis, no signs of instability are present for 9 above 
a certain threshold angle (about 9 ~ 0.23 deg in Fig. 
0 , but a definitive response to the question of the ab¬ 
solute stabilization of NL-UBBs cannot be given. The 
trend of — ImK with increasing cone angle suggests an 
exponential decay. With increasing intensity, Kerr non¬ 
linearity renders NL-UBBs increasingly unstable at first, 
but the increasing NLLs has an opposite stabilizing effect 
at higher intensities. Stabilization by NLL appears to be 
complete above a certain cone angle (about 0 ~ 0.15 deg 
in Fig. 0. It is interesting that below this angle, once 
the Kerr-induced unstable mode disappears by the ac¬ 
tion of NLLs (at about Iq = 18 TW/cm^), an underlying 
unstable mode with a different oscillation frequency Ken 
becomes dominant. 

In connection with Figs. 0a) and 0a), the steady 
propagation regime after the axicon can be seen to be 
associated with the absence of unstable modes of the at¬ 
tracting NL-UBB with 9 = 0.15 and Iq = 28 TW/cm^ 
that tends to be formed at the center of the Bessel zone. 
Even if repeated self-focusing cycles are observed before 
the focus of the axicon for these small cone angles (cy¬ 


cles that may be much more pronounced), and before 
the waist of the Bessel-Gauss beam, the input radiation 
is pushed stably towards the NL-UBB about the center of 
the Bessel zone. In Fig. 0c) and Fig. 0b), instead, the 
(weak) unstable regime appears to reflect the instability 
of the attracting NL-UBB with 9 = 0.15 and Iq = 12 
TW/cm^. The gain is indeed low [Fig. 0c)] (compared 
to next situations below), and the oscillation frequency 
about the focus or waist is seen to coincide with the os¬ 
cillation frequency Re/t of the dominant unstable mode 
of the attracting NL-UBB [Fig. 0d)]. 

The connection between unsteady Bessel propaga¬ 
tion regime and instability of the attracting NL-UBB is 
clearer in situations of higher gain. For two NL-UBBs 
with increasing gain. Fig. 0a) and (b) shows the growth 
of the respective dominant unstable modes with propa¬ 
gation distance. These small modes develop into large, 
periodic (but no longer harmonic) perturbation regimes, 
that gradually turn into quasi-periodic, and eventually 
into chaos. In all cases we have studied, this process is 
found to be faster as the gain —ImK triggering this pro¬ 
cess is higher. Also, the oscillation frequency in the large, 
periodic perturbation regime is close to but slightly lower 
than ReK. In simulations as those of Fig. 0 NL-UBBs 
are directly launched into the medium and the dominant, 
unstable modes of each NL-UBB are found to emerge 
spontaneously from numerical noise with the gain and os¬ 
cillation frequency predicted by the linearized instability 
analysis [left part in Figs. 0a) and (b)]. To simulate the 



































































































propagation of ideal, non-truncated NL-UBBs, we use 
the procedure of replacing the propagated field at each 
axial numerical step of propagation with the initial NL- 
UBB in a narrow annulus touching the end of the (quite 
large) numerical radial grid. This procedure is justified 
since no dynamics is expected to take place in the linear 
asymptotic tails. In Figs. [HJc-f), phase spaces I-dl/dz 
(/ on-axis intensity) in relevant propagation intervals are 
shown. In the case of lower gain, the large perturbation 
regime remains periodic for a considerable propagation 
distance [Fig. HDc)], whereas in the case of higher gain, it 
becomes quasi-periodic from the beginning of the large 
perturbation regime [Fig. EDd)] and enters sooner into 
chaos. The phase spaces up to the longest propagation 
distance [Figs. |6l[e) and (f)] evidence that NL-UBBs are 
actually chaotic attractors, whose morphology depends 
on the specific NL-UBB. 

On the other hand, Figs. |71[a) and (b) show the on-axis 
intensities after an axicon illuminated with two Gaussian 
beams of width and intensities such that the attractors 
are the two NL-UBBs analyzed above. In the case of Fig. 
|71[a) with lower gain, the oscillation frequency about the 
focus of the axicon coincides with that of the large, pe¬ 
riodic, perturbation regime. Indeed the structure of the 
oscillations in the phase space of Fig. me) about the fo¬ 
cus of the axicon mimics the structure of the anharmonic 
oscillations in the periodic perturbation regime of the at¬ 
tracting NL-UBB in Fig. ETc). Small differences orig¬ 
inate from the slow decay of intensity along the Bessel 
zone of the finite-power BB. The structure of the phase 
space in the whole Bessel zone (not shown) does not re¬ 
produce the morphology of the chaotic attractor in Fig. 
Ele). In the case of Fig. |71[b) with higher gain, the on- 
axis intensity in a considerable part of the Bessel zone 
exhibits a highly disordered dynamics. Comparison of 
the morphology of the phase space in Fig. |71[d) for the 
whole Bessel zone with that of the attracting NL-UBB 
in Fig. EKf) evidences that the dynamics after the axi¬ 
con is reproducing the chaotic dynamics about the ideal 
NL-UBB chaotic attractor. 


VI. CONCLUSIONS 

From a series of diagnostic numerical simulations, we 
have extracted the underlying laws governing the spatial 
dynamics of the light beam emerging from an axicon, 
and entering a medium where self-focusing Kerr effect 
and multiphoton absorption are relevant. If as pointed 


out, temporal and plasma effects play a secondary role in 
determining the spatial dynamics in filamentation with 
BBs, these laws provide an unified understanding of the 
different Bessel filamentation regimes described previ¬ 
ously. In a few words, the nonlinear propagation is deter¬ 
mined by an attracting NL-UBB and its stability proper¬ 
ties under small perturbations. The attracting NL-UBB 
is that whose inward Hankel amplitude equals the am¬ 
plitude of the BB that the BB generator would create 
at the center of the Bessel zone in linear propagation. 
We have derived an approximate analytical expression 
that determines the attracting NL-UBB given the opti¬ 
cal properties of the medium, the cone angle, and the 
intensity of the linear BB (or equivalently, the axicon 
base angle and the input Gaussian width and intensity). 
Steady/unsteady propagation regimes are shown to cor¬ 
respond to stability/instability of the attracting NL-UBB 
under small radial perturbations, i. e., to the existence of 
a small unstable radial mode that tends to grow exponen¬ 
tially. We have performed an extensive stability analysis 
under small radial perturbations that put in quantitative 
terms the stabilization effect of increasing the cone an¬ 
gle and the intensity. In case of instability under small 
perturbations, NL-UBBs are seen to develop a large per¬ 
turbation and chaotic regimes with increasing propaga¬ 
tion distances. In the Bessel zone after the axicon, and 
depending on how large the gain of the small unstable 
mode of the attracting NL-UBB is, the unsteady dynam¬ 
ics reproduces the dynamics of the small perturbation, 
large, or chaotic perturbation regimes of the attracting 
NL-UBB. Though a direct relation with increasing gain 
is obvious, further research would be needed to specify 
in more quantitative terms the particular perturbation 
regime (small, large or chaotic) of the attracting NL-UBB 
that is observed in the Bessel zone. We have restricted 
ourselves to vortex-less NL-UBBs, but the generality of 
these ideas suggests a relatively simple generalization to 
axicon-generated vortex NL-UBB, in which case not only 
radial instability but also azimuthal instability should be 
taken into account. 
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